rm(list=ls())
library(foreign)
library(MASS)
library(survival)

#Load Dataset
dat <- read.dta("use.dta")

#Re-name counting process variables to R format
dat$t0 <- dat$"_t0"
dat$t <- dat$"_t"

#Generate log-time interactions
dat$eth_tlog <- dat$eth*log(dat$t)
dat$rel_tlog <- dat$rel*log(dat$t)
dat$mtn_tlog <- dat$mtn*log(dat$t)
dat$pol_tlog <- dat$polity2*log(dat$t)

###################
##### Table 2 #####
###################

gv <- coxph(Surv(t0, t, govvictory) ~ indirebel_dec + indigov_dec + bruterebel_dec + brutegov_dec + brutegov_dec_tlog + coercrebel_dec + coercgov_dec + oppose_dec + rel + intens + mtn + eth + coldwar + polity2 + oppose_dec + num_dyad + opp_lt + cluster(conflict), data=dat, method="efron")
summary(gv)

rv <- coxph(Surv(t0, t, rebvictory) ~ indirebel_dec + indigov_dec + bruterebel_dec + bruterebel_dec_lt + brutegov_dec + coercrebel_dec + coercgov_dec + rel + intens + mtn + eth + coldwar + polity2 + oppose_dec + num_dyad + opp_lt + eth_tlog + cluster(conflict), data=dat, method="efron")
summary(rv)

ns <- coxph(Surv(t0, t, negsett) ~ indirebel_dec + indirebel_dec_lt + indigov_dec + bruterebel_dec + brutegov_dec + brutegov_dec_tlog + coercrebel_dec + rel + intens + mtn + eth + coldwar + polity2 + oppose_dec + num_dyad + opp_lt + rel_tlog + mtn_tlog + pol_tlog + eth_tlog + cluster(conflict), data=dat, method="efron")
summary(ns)


############################
######### Figure 1 #########
############################

###Simluating Coefficient for Indirect Rebel*Log Time on Negotiated Settlement
sims<-10000
set.seed(031415)
sim.betas<-mvrnorm(sims,mu=ns$coef,Sigma= ns$var)
b1<-sim.betas[,1]
b2<-sim.betas[,2]

t.sim<-seq(1, 350, by=1)
out<-matrix(0, nrow=length(t.sim), ncol=3)

for(i in 1:length(t.sim)){
	simfit<- ((b1+ (b2*log(t.sim[i]))))
	out[i,]<-c(mean(simfit), quantile(simfit, c(.025, .975)))
	}

#Find where the effect loses significance - loses significance at 44 months
cbind(t.sim[39:49],out[39:49,2:3])

#Plotting the Simulation
plot(t.sim, out[,1], type="n", xlab=list("Time Since Conflict Onset (months)",cex=.9), ylab=list("Combined Coefficient", cex=.9), main="Negotiated Settlement", xlim=c(0,250))
lines(t.sim, out[,1], lty=1)
lines(t.sim, out[,2], lty=2)
lines(t.sim, out[,3], lty=2)
abline(h=0, lwd=3)
indireb_t <- dat$indirebel*dat$t
rug(indireb_t, lwd=1)

#Draw line where lower confidence bound crosses 0
abline(v=44, col="grey")


############################
######### Figure 2 #########
############################

###Simluating Coefficient for Direct-Conventional Rebel*Log Time on Rebel Victory
sims<-10000
set.seed(031415)
sim.betas<-mvrnorm(sims,mu=rv$coef,Sigma=rv$var)
b2<-sim.betas[,3]
b3<-sim.betas[,4]

t.sim<-seq(1, 350, by=1)
out<-matrix(0, nrow=length(t.sim), ncol=3)

for(i in 1:length(t.sim)){
	simfit<- ((b2+ (b3*log(t.sim[i]))))
	out[i,]<-c(mean(simfit), quantile(simfit, c(.025, .975)))
	}

#Find where the effect loses statistical significance - loses significance at 25 months
cbind(t.sim[20:30], out[20:30,2:3])

#Plotting the Simulation
plot(t.sim, out[,1], type="n", xlab=list("Time Since Conflict Onset (months)", cex=1), ylab=list("Combined Coefficient", cex=1), main="Rebel Victory", xlim=c(0,250))
lines(t.sim, out[,1], lty=1)
lines(t.sim, out[,2], lty=2)
lines(t.sim, out[,3], lty=2)
abline(h=0, lwd=3)
bruterebel_t <- dat$bruterebel * dat$t
rug(bruterebel_t)

#Draw line where lower confidence bound crosses 0
abline(v=25, col="grey")


#############################
######### Figure 3a #########
#############################

###Simluating Coefficient for Direct-Conventional Government*Log Time on Government Victory
sims<-10000
set.seed(031415)
sim.betas<-mvrnorm(sims,mu=gv$coef,Sigma=gv$var)
b1<-sim.betas[,4]
b2<-sim.betas[,5]

t.sim<-seq(1, 350, by=1)
out<-matrix(0, nrow=length(t.sim), ncol=3)

for(i in 1:length(t.sim)){
	simfit<- ((b1+ (b2*log(t.sim[i]))))
	out[i,]<-c(mean(simfit), quantile(simfit, c(.025, .975)))
	}

#Find where the effect loses significance - loses significance at 10 months
cbind(t.sim[5:15],out[5:15,2:3])

#Find where the effect re-gains significance - re-gains significance at 111 months
cbind(t.sim[106:116],out[106:116,2:3])


#Plotting the Simulation
plot(t.sim, out[,1], type="n", xlab=list("Time Since Conflict Onset (months)", cex=1), ylab=list("Combined Coefficient", cex=1), main= list("Government Victory", cex=1), ylim=c(-10,10), xlim=c(0,250))
lines(t.sim, out[,1], lty=1)
lines(t.sim, out[,2], lty=2)
lines(t.sim, out[,3], lty=2)
abline(h=0, lwd=3)
brutegov_t <- dat$brutegov*dat$t
rug(brutegov_t, lwd=1)
#Draw lines where the effect gains and loses statistical significance
abline(v=10, col="grey")
abline(v=111, col="grey")


#############################
######### Figure 3b #########
#############################

###Simluating Coefficient for Direct-Conventional Government*Log Time on Negotiated Settlement
sims<-10000
set.seed(031415)
sim.betas<-mvrnorm(sims,mu=ns$coef,Sigma=ns$var)
b1<-sim.betas[,5]
b2<-sim.betas[,6]

t.sim<-seq(1, 350, by=1)
out<-matrix(0, nrow=length(t.sim), ncol=3)

for(i in 1:length(t.sim)){
	simfit<- ((b1+ (b2*log(t.sim[i]))))
	out[i,]<-c(mean(simfit), quantile(simfit, c(.025, .975)))
	}

#Find where the effect gains significance - gains significance at 18 months
cbind(t.sim[13:24],out[13:24,2:3])

#The effect is briefly significant at t=1, as noted in footnote 26
cbind(t.sim[0:5],out[0:5,2:3])

#Plotting the Simulation
plot(t.sim, out[,1], type="n", xlab=list("Time Since Conflict Onset (months)",cex=.9), ylab=list("Combined Coefficient", cex=.9), main="Negotiated Settlement", xlim=c(0,250))
lines(t.sim, out[,1], lty=1)
lines(t.sim, out[,2], lty=2)
lines(t.sim, out[,3], lty=2)
abline(h=0, lwd=3)
brutegov_t <- dat$brutegov*dat$t
rug(brutegov_t, lwd=1)
#Draw lines where the effect gains statistical significance
abline(v=18, col="grey")

